library(data.table)
library(readr)
library(dplyr)
library(knitr)
library(kableExtra)
library(Seurat)
#seurobj <- readRDS('output/seurat_objects/180831/10x-180831')
gene_annotations <- fread('../tables/tables_paper/all_tables/genes_biomart.txt', sep='\t', quote="", header=T)

rr ids2symbols <- read.table(‘tables/tables_paper/all_tables/10x-180831-geneids-genesymbols.tsv’, header=T)

BEAM clusters

combined <- list()

for (i in 1:6){
  beam_genelist <- read.table(paste('tables/tables_paper/all_tables/BEAM/heatmap_logFC0.3_ncluster6/genelist_cluster', i, '.tsv', sep=''), sep='\t', header=T)
  beam_genelist['ensembl_gene_id'] <- ids2symbols$ensembl_gene_id[match(beam_genelist$gene_short_name, ids2symbols$gene_symbol)]
  beam_genelist['cluster'] <- i
  
  #add geneinfo
  merged <- merge(beam_genelist, gene_annotations[,c('Gene stable ID', 'Gene type', 'GO term accession', 'GO term name')], by.y='Gene stable ID', by.x='ensembl_gene_id')
  
  #filter for tf's and gene type != protein coding
  filtered <- merged[c(which(c(merged$`GO term name`) %like% 'transcription factor'),
                     which(merged$`GO term definition` %like% 'transcription factor'),
                     which(merged$`Gene type` != 'protein_coding')),]
  
  #filter out duplicate genes
  filtered <- filtered[!duplicated(filtered$gene_short_name),]
  
  combined[[i]] <- filtered
}

combined <- do.call("rbind", combined)
kable(combined[c('gene_short_name', 'qval', 'avgLogFC_State2_State3', 'cluster', 'Gene type', 'GO term accession', 'GO term name')]) %>%
  kable_styling(bootstrap_options = "striped")
<table class="table table-striped" style="margin-left: auto; margin-right: auto;">
 <thead>
  <tr>
   <th style="text-align:left;">   </th>
   <th style="text-align:left;"> gene_short_name </th>
   <th style="text-align:right;"> qval </th>
   <th style="text-align:right;"> avgLogFC_State2_State3 </th>
   <th style="text-align:right;"> cluster </th>
   <th style="text-align:left;"> Gene type </th>
   <th style="text-align:left;"> GO term accession </th>
   <th style="text-align:left;"> GO term name </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;"> 1614 </td>
   <td style="text-align:left;"> C1QBP </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.4210100 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0008134 </td>
   <td style="text-align:left;"> transcription factor binding </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4113 </td>
   <td style="text-align:left;"> YWHAZ </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3358994 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0008134 </td>
   <td style="text-align:left;"> transcription factor binding </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4589 </td>
   <td style="text-align:left;"> PHB </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3788402 </td>
   <td style="text-align:right;"> 1 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 238 </td>
   <td style="text-align:left;"> MLXIPL </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3673746 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 754 </td>
   <td style="text-align:left;"> NR1H3 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3569965 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0003700 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1808 </td>
   <td style="text-align:left;"> SREBF1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.4302819 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4371 </td>
   <td style="text-align:left;"> CAT </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3321904 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0032088 </td>
   <td style="text-align:left;"> negative regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 5709 </td>
   <td style="text-align:left;"> PPARG </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3226131 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0001103 </td>
   <td style="text-align:left;"> RNA polymerase II repressing transcription factor binding </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 7406 </td>
   <td style="text-align:left;"> CD36 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 1.4585314 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051092 </td>
   <td style="text-align:left;"> positive regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 11046 </td>
   <td style="text-align:left;"> MESP1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.4067898 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0003700 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 11896 </td>
   <td style="text-align:left;"> FZD4 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3027719 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051091 </td>
   <td style="text-align:left;"> positive regulation of DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 12856 </td>
   <td style="text-align:left;"> PTMA </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.6008086 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0033613 </td>
   <td style="text-align:left;"> activating transcription factor binding </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 14113 </td>
   <td style="text-align:left;"> CEBPA </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3466449 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 13969 </td>
   <td style="text-align:left;"> EPHA1-AS1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3000552 </td>
   <td style="text-align:right;"> 2 </td>
   <td style="text-align:left;"> antisense </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 264 </td>
   <td style="text-align:left;"> SNAI2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5313284 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1320 </td>
   <td style="text-align:left;"> EGR1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4104048 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1409 </td>
   <td style="text-align:left;"> CLU </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.7825588 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051092 </td>
   <td style="text-align:left;"> positive regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2373 </td>
   <td style="text-align:left;"> ZEB1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4307388 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2683 </td>
   <td style="text-align:left;"> ZFP36L2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3991944 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0003700 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2993 </td>
   <td style="text-align:left;"> PPAP2B </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5215079 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051091 </td>
   <td style="text-align:left;"> positive regulation of DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3109 </td>
   <td style="text-align:left;"> OSR2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5893658 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3745 </td>
   <td style="text-align:left;"> FOS </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.6405352 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051090 </td>
   <td style="text-align:left;"> regulation of DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3982 </td>
   <td style="text-align:left;"> CD34 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4275094 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0008134 </td>
   <td style="text-align:left;"> transcription factor binding </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4254 </td>
   <td style="text-align:left;"> ZFP36L1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5867891 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0003700 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4537 </td>
   <td style="text-align:left;"> MIR24-2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3113038 </td>
   <td style="text-align:right;"> 3 </td>
   <td style="text-align:left;"> antisense </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2073 </td>
   <td style="text-align:left;"> PTGIS </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3420740 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0032088 </td>
   <td style="text-align:left;"> negative regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2661 </td>
   <td style="text-align:left;"> CYP1B1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.6809188 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0032088 </td>
   <td style="text-align:left;"> negative regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3478 </td>
   <td style="text-align:left;"> ARID5B </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4524769 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3580 </td>
   <td style="text-align:left;"> CARHSP1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3443580 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 5930 </td>
   <td style="text-align:left;"> TMSB4X </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3553073 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0032088 </td>
   <td style="text-align:left;"> negative regulation of NF-kappaB transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 4861 </td>
   <td style="text-align:left;"> MIR4435-1HG </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4703532 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> lincRNA </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 5050 </td>
   <td style="text-align:left;"> ZFAS1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3843600 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> antisense </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 6029 </td>
   <td style="text-align:left;"> LINC00152 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5044662 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> lincRNA </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 6067 </td>
   <td style="text-align:left;"> RP11-14N7.2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.3392808 </td>
   <td style="text-align:right;"> 4 </td>
   <td style="text-align:left;"> lincRNA </td>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;">  </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 411 </td>
   <td style="text-align:left;"> PRRX1 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.5139881 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1868 </td>
   <td style="text-align:left;"> DDR2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4279312 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0051091 </td>
   <td style="text-align:left;"> positive regulation of DNA-binding transcription factor activity </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2060 </td>
   <td style="text-align:left;"> JUNB </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.7652085 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 2452 </td>
   <td style="text-align:left;"> TCF4 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4466219 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0005667 </td>
   <td style="text-align:left;"> transcription factor complex </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 3593 </td>
   <td style="text-align:left;"> VGLL3 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> -0.4709351 </td>
   <td style="text-align:right;"> 5 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1562 </td>
   <td style="text-align:left;"> BOLA3 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.3568609 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0000981 </td>
   <td style="text-align:left;"> DNA-binding transcription factor activity, RNA polymerase II-specific </td>
  </tr>
  <tr>
   <td style="text-align:left;"> 1608 </td>
   <td style="text-align:left;"> PRDX2 </td>
   <td style="text-align:right;"> 0 </td>
   <td style="text-align:right;"> 0.5220889 </td>
   <td style="text-align:right;"> 6 </td>
   <td style="text-align:left;"> protein_coding </td>
   <td style="text-align:left;"> GO:0032088 </td>
   <td style="text-align:left;"> negative regulation of NF-kappaB transcription factor activity </td>
  </tr>
</tbody>
</table>

rr symbols2ids <- hs_add_gene_symbol_from_ensembl_ids(data.frame(colname_geneids_from=rownames(seurobj@data), colname_geneids_to=NA))

Error in `[.data.frame`(df, , colname_geneids_from) : 
  undefined columns selected
plots_upper <- FeaturePlot(seurobj, features.plot=as.vector(combined_upper_branch$gene_short_name), nCol=4, cols.use=c('gray', 'blue'), no.legend=F, no.axes=T, do.return=T)

plots_upper_edited <- list()
for (p in names(plots_upper)){
  plots_upper_edited[[p]] <- plots_upper[[p]] + scale_color_gradient(name='Expr.', low='gray', high='blue', guide='colorbar') + theme(legend.title=element_text(size=9), legend.text=element_text(size=9), legend.key.height = unit(0.6, 'cm'), legend.key.width=unit(0.2, 'cm'))
}
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
grid_upper <- plot_grid(plotlist=plots_upper_edited, ncol=4)
grid_upper

#save_plot('figures/figures_paper/supplementary_figures/10_tf-analysis/TFs_upper_branch.pdf', grid_upper, base_height=9.14, base_width=12)
#save_plot('figures/figures_paper/supplementary_figures/10_tf-analysis/TFs_upper_branch.png', grid_upper, base_height=9.14, base_width=12)
plots_lower <- FeaturePlot(seurobj, features.plot=as.vector(combined_lower_branch$gene_short_name), nCol=3, cols.use=c('gray', 'blue'), no.legend=F, no.axes=T, do.return=T)

plots_lower_edited <- list()
for (p in names(plots_lower)){
  plots_lower_edited[[p]] <- plots_lower[[p]] + scale_color_gradient(name='Expr.', low='gray', high='blue', guide='colorbar') + theme(legend.title=element_text(size=9), legend.text=element_text(size=9), legend.key.height = unit(0.6, 'cm'), legend.key.width=unit(0.2, 'cm'))
}
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
Scale for 'colour' is already present. Adding another scale for 'colour',
which will replace the existing scale.
grid_lower <- plot_grid(plotlist=plots_lower_edited, ncol=4)

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCgpgYGB7ciBtZXNzYWdlPUZ9CmxpYnJhcnkoZGF0YS50YWJsZSkKbGlicmFyeShyZWFkcikKbGlicmFyeShkcGx5cikKbGlicmFyeShrbml0cikKbGlicmFyeShrYWJsZUV4dHJhKQpsaWJyYXJ5KFNldXJhdCkKCnNldXJvYmogPC0gcmVhZFJEUygnb3V0cHV0L3NldXJhdF9vYmplY3RzLzE4MDgzMS8xMHgtMTgwODMxJykKZ2VuZV9hbm5vdGF0aW9ucyA8LSBmcmVhZCgndGFibGVzL3RhYmxlc19wYXBlci9hbGxfdGFibGVzL2dlbmVzX2Jpb21hcnQudHh0Jywgc2VwPSdcdCcsIHF1b3RlPSIiLCBoZWFkZXI9VCkKYGBgCgoKYGBge3J9CmlkczJzeW1ib2xzIDwtIHJlYWQudGFibGUoJ3RhYmxlcy90YWJsZXNfcGFwZXIvYWxsX3RhYmxlcy8xMHgtMTgwODMxLWdlbmVpZHMtZ2VuZXN5bWJvbHMudHN2JywgaGVhZGVyPVQpCmBgYAoKCiNCRUFNIGNsdXN0ZXJzCgpgYGB7cn0KY29tYmluZWQgPC0gbGlzdCgpCgpmb3IgKGkgaW4gMTo2KXsKICBiZWFtX2dlbmVsaXN0IDwtIHJlYWQudGFibGUocGFzdGUoJ3RhYmxlcy90YWJsZXNfcGFwZXIvYWxsX3RhYmxlcy9CRUFNL2hlYXRtYXBfbG9nRkMwLjNfbmNsdXN0ZXI2L2dlbmVsaXN0X2NsdXN0ZXInLCBpLCAnLnRzdicsIHNlcD0nJyksIHNlcD0nXHQnLCBoZWFkZXI9VCkKICBiZWFtX2dlbmVsaXN0WydlbnNlbWJsX2dlbmVfaWQnXSA8LSBpZHMyc3ltYm9scyRlbnNlbWJsX2dlbmVfaWRbbWF0Y2goYmVhbV9nZW5lbGlzdCRnZW5lX3Nob3J0X25hbWUsIGlkczJzeW1ib2xzJGdlbmVfc3ltYm9sKV0KICBiZWFtX2dlbmVsaXN0WydjbHVzdGVyJ10gPC0gaQogIAogICNhZGQgZ2VuZWluZm8KICBtZXJnZWQgPC0gbWVyZ2UoYmVhbV9nZW5lbGlzdCwgZ2VuZV9hbm5vdGF0aW9uc1ssYygnR2VuZSBzdGFibGUgSUQnLCAnR2VuZSB0eXBlJywgJ0dPIHRlcm0gYWNjZXNzaW9uJywgJ0dPIHRlcm0gbmFtZScpXSwgYnkueT0nR2VuZSBzdGFibGUgSUQnLCBieS54PSdlbnNlbWJsX2dlbmVfaWQnKQogIAogICNmaWx0ZXIgZm9yIHRmJ3MgYW5kIGdlbmUgdHlwZSAhPSBwcm90ZWluIGNvZGluZwogIGZpbHRlcmVkIDwtIG1lcmdlZFtjKHdoaWNoKGMobWVyZ2VkJGBHTyB0ZXJtIG5hbWVgKSAlbGlrZSUgJ3RyYW5zY3JpcHRpb24gZmFjdG9yJyksCiAgICAgICAgICAgICAgICAgICAgIHdoaWNoKG1lcmdlZCRgR08gdGVybSBkZWZpbml0aW9uYCAlbGlrZSUgJ3RyYW5zY3JpcHRpb24gZmFjdG9yJyksCiAgICAgICAgICAgICAgICAgICAgIHdoaWNoKG1lcmdlZCRgR2VuZSB0eXBlYCAhPSAncHJvdGVpbl9jb2RpbmcnKSksXQogIAogICNmaWx0ZXIgb3V0IGR1cGxpY2F0ZSBnZW5lcwogIGZpbHRlcmVkIDwtIGZpbHRlcmVkWyFkdXBsaWNhdGVkKGZpbHRlcmVkJGdlbmVfc2hvcnRfbmFtZSksXQogIAogIGNvbWJpbmVkW1tpXV0gPC0gZmlsdGVyZWQKfQoKY29tYmluZWQgPC0gZG8uY2FsbCgicmJpbmQiLCBjb21iaW5lZCkKYGBgCgpgYGB7cn0Ka2FibGUoY29tYmluZWRbYygnZ2VuZV9zaG9ydF9uYW1lJywgJ3F2YWwnLCAnYXZnTG9nRkNfU3RhdGUyX1N0YXRlMycsICdjbHVzdGVyJywgJ0dlbmUgdHlwZScsICdHTyB0ZXJtIGFjY2Vzc2lvbicsICdHTyB0ZXJtIG5hbWUnKV0pICU+JQogIGthYmxlX3N0eWxpbmcoYm9vdHN0cmFwX29wdGlvbnMgPSAic3RyaXBlZCIpCmBgYAoKCmBgYHtyfQojd3JpdGUudGFibGUoY29tYmluZWQsIGZpbGU9Jy4uL3RhYmxlcy90YWJsZXNfcGFwZXIvYWxsX3RhYmxlcy9CRUFNL2hlYXRtYXBfbG9nRkMwLjNfbmNsdXN0ZXI2L2dlbmVsaXN0X1RGcy50c3YnLCBzZXA9J1x0JywgcXVvdGU9Riwgcm93Lm5hbWVzPUYpCiN3cml0ZS50YWJsZShjb21iaW5lZCwgZmlsZT0nLi4vdGFibGVzL3RhYmxlc19wYXBlci9zdXBwbGVtZW50YXJ5X3RhYmxlcy9CRUFNX2xvZ0ZDMC4zX25jbHVzdGVyczYvZ2VuZWxpc3RfVEZzLnRzdicsIHNlcD0nXHQnLCBxdW90ZT1GLCByb3cubmFtZXM9RikKYGBgCgpgYGB7cn0KY29tYmluZWRfdXBwZXJfYnJhbmNoIDwtIGNvbWJpbmVkW2NvbWJpbmVkJGF2Z0xvZ0ZDX1N0YXRlMl9TdGF0ZTMgPiAwLF0KY29tYmluZWRfbG93ZXJfYnJhbmNoIDwtIGNvbWJpbmVkW2NvbWJpbmVkJGF2Z0xvZ0ZDX1N0YXRlMl9TdGF0ZTMgPCAwLF0KYGBgCgpgYGB7ciBmaWcuc2hvdz0naGlkZScsIG1lc3NhZ2U9Rn0KcGxvdHNfdXBwZXIgPC0gRmVhdHVyZVBsb3Qoc2V1cm9iaiwgZmVhdHVyZXMucGxvdD1hcy52ZWN0b3IoY29tYmluZWRfdXBwZXJfYnJhbmNoJGdlbmVfc2hvcnRfbmFtZSksIG5Db2w9NCwgY29scy51c2U9YygnZ3JheScsICdibHVlJyksIG5vLmxlZ2VuZD1GLCBuby5heGVzPVQsIGRvLnJldHVybj1UKQoKcGxvdHNfdXBwZXJfZWRpdGVkIDwtIGxpc3QoKQoKZm9yIChwIGluIG5hbWVzKHBsb3RzX3VwcGVyKSl7CiAgcGxvdHNfdXBwZXJfZWRpdGVkW1twXV0gPC0gcGxvdHNfdXBwZXJbW3BdXSArIHNjYWxlX2NvbG9yX2dyYWRpZW50KG5hbWU9J0V4cHIuJywgbG93PSdncmF5JywgaGlnaD0nYmx1ZScsIGd1aWRlPSdjb2xvcmJhcicpICsgdGhlbWUobGVnZW5kLnRpdGxlPWVsZW1lbnRfdGV4dChzaXplPTkpLCBsZWdlbmQudGV4dD1lbGVtZW50X3RleHQoc2l6ZT05KSwgbGVnZW5kLmtleS5oZWlnaHQgPSB1bml0KDAuNiwgJ2NtJyksIGxlZ2VuZC5rZXkud2lkdGg9dW5pdCgwLjIsICdjbScpKQp9CmdyaWRfdXBwZXIgPC0gcGxvdF9ncmlkKHBsb3RsaXN0PXBsb3RzX3VwcGVyX2VkaXRlZCwgbmNvbD00KQpgYGAKCgpgYGB7ciBmaWcuaGVpZ2h0ID0gMTAsIGZpZy53aWR0aCA9IDEyLCBmaWcuYWxpZ24gPSAiY2VudGVyIn0KZ3JpZF91cHBlcgojc2F2ZV9wbG90KCdmaWd1cmVzL2ZpZ3VyZXNfcGFwZXIvc3VwcGxlbWVudGFyeV9maWd1cmVzLzEwX3RmLWFuYWx5c2lzL1RGc191cHBlcl9icmFuY2gucGRmJywgZ3JpZF91cHBlciwgYmFzZV9oZWlnaHQ9OS4xNCwgYmFzZV93aWR0aD0xMikKI3NhdmVfcGxvdCgnZmlndXJlcy9maWd1cmVzX3BhcGVyL3N1cHBsZW1lbnRhcnlfZmlndXJlcy8xMF90Zi1hbmFseXNpcy9URnNfdXBwZXJfYnJhbmNoLnBuZycsIGdyaWRfdXBwZXIsIGJhc2VfaGVpZ2h0PTkuMTQsIGJhc2Vfd2lkdGg9MTIpCmBgYAoKYGBge3IgZmlnLnNob3c9J2hpZGUnLCBtZXNzYWdlPUZ9CnBsb3RzX2xvd2VyIDwtIEZlYXR1cmVQbG90KHNldXJvYmosIGZlYXR1cmVzLnBsb3Q9YXMudmVjdG9yKGNvbWJpbmVkX2xvd2VyX2JyYW5jaCRnZW5lX3Nob3J0X25hbWUpLCBuQ29sPTMsIGNvbHMudXNlPWMoJ2dyYXknLCAnYmx1ZScpLCBuby5sZWdlbmQ9Riwgbm8uYXhlcz1ULCBkby5yZXR1cm49VCkKCnBsb3RzX2xvd2VyX2VkaXRlZCA8LSBsaXN0KCkKCmZvciAocCBpbiBuYW1lcyhwbG90c19sb3dlcikpewogIHBsb3RzX2xvd2VyX2VkaXRlZFtbcF1dIDwtIHBsb3RzX2xvd2VyW1twXV0gKyBzY2FsZV9jb2xvcl9ncmFkaWVudChuYW1lPSdFeHByLicsIGxvdz0nZ3JheScsIGhpZ2g9J2JsdWUnLCBndWlkZT0nY29sb3JiYXInKSArIHRoZW1lKGxlZ2VuZC50aXRsZT1lbGVtZW50X3RleHQoc2l6ZT05KSwgbGVnZW5kLnRleHQ9ZWxlbWVudF90ZXh0KHNpemU9OSksIGxlZ2VuZC5rZXkuaGVpZ2h0ID0gdW5pdCgwLjYsICdjbScpLCBsZWdlbmQua2V5LndpZHRoPXVuaXQoMC4yLCAnY20nKSkKfQpncmlkX2xvd2VyIDwtIHBsb3RfZ3JpZChwbG90bGlzdD1wbG90c19sb3dlcl9lZGl0ZWQsIG5jb2w9NCkKYGBgCgpgYGB7ciwgZmlnLmhlaWdodCA9IDE1LCBmaWcud2lkdGggPSAxMiwgZmlnLmFsaWduID0gImNlbnRlciJ9CmdyaWRfbG93ZXIKI3NhdmVfcGxvdCgnZmlndXJlcy9maWd1cmVzX3BhcGVyL3N1cHBsZW1lbnRhcnlfZmlndXJlcy8xMF90Zi1hbmFseXNpcy9URnNfbG93ZXJfYnJhbmNoLnBkZicsIGdyaWRfbG93ZXIsIGJhc2VfaGVpZ2h0PTE2LCBiYXNlX3dpZHRoPTEyKQojc2F2ZV9wbG90KCdmaWd1cmVzL2ZpZ3VyZXNfcGFwZXIvc3VwcGxlbWVudGFyeV9maWd1cmVzLzEwX3RmLWFuYWx5c2lzL1RGc19sb3dlcl9icmFuY2gucG5nJywgZ3JpZF9sb3dlciwgYmFzZV9oZWlnaHQ9MTYsIGJhc2Vfd2lkdGg9MTIpCmBgYAoK